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ABSTRACT 

This  report  describes  work  carried  out  under  contract  SPC-96-4054  ‘‘Optimal  Shack- 
Hartmann  wavefront  sensing  for  low  light  levels” .  The  aim  of  this  work  has  been  to  identify, 
develop  and,  where  appropriate,  implement  methods  for  improving  wavefront  estimates  at 
low  light  levels.  The  report  is  divided  into  three  main  sections  followed  by  a  summary  and 
appendices.  Each  section  addresses  a  major  aspect  of  the  work  car  lied  out. 

•  The  first  section  considers  the  question  of  how  to  optimally  choose  the  size  of  Hartmann 
subapertures  for  a  given  signal  level  so  as  to  achieve  the  optimal  wavefront  estimate.  After 
a  short  introduction,  the  main  body  of  this  section  is  given  in  the  form  of  a  paper  recently 
submitted  for  publication  in  Optics  Communications  entitled  "Optimal  Hartmann  sensing 
at  low  fight  levels” . 

•  Section  2  describes  the  use  of  a  Gram-Charfier  matched  filter  to  yield  improved  sub- 
aperture  centroid  estimates.  The  importance  of  the  centroiding  operation  in  wavefront 
sensing  and  adaptive  optics  applications  is  first  described  and  the  theory  and  results  ob¬ 
tained  with  this  matched  filter  are  presented  in  the  form  of  a  paper  recently  submitted  and 
accepted  by  Optics  Letters  -  “Gram-Charfier  matched  filter  for  Shack- Hartmann  sensing 
at  low  light  levels” . 

•  The  third  section  discusses  the  use  of  both  Bayesian  estimators  and  modal  projection 
estimation  (combined  with  unconventional  Hartmann  sensor  geometry)  as  a  means  to 
reduce  noise  propagation  in  wavefront  sensing  applications  operating  at  low  fight  level. 
Theory  is  presented  together  with  computer  simulated  results  for  both  techniques. 

•  The  results  achieved  and  understanding  gained  are  given  in  the  summary. 


19971209  025 


1.  OPTIMAL  SUB- APERTURE  SIZE  FOR  HARTMANN  SENSORS 

When  Hartmann  sensors  are  required  to  operate  in  low  signal-noise  environments, 
there  is  an  inherent  trade-off  regarding  how  the  available  signal  photons  should  be  dis¬ 
tributed”.  The  Hartmann  sensor  is  designed  to  provide  estimates  of  the  local  slope  of  a 
wavefront  averaged  over  a  number  of  sample  areas  in  the  pupil  plane.  From  these  averaged 
slope  measurements,  the  wavefront  estimate  must  directly  (or  indirectly)  be  calculated. 
The  linear  estimation  problem  that  arises  has  been  derived  and  discussed  in  section  1  of 
the  interim  report1.  The  key  point  here  is  that  whilst  choice  of  a  large  number  of  sub¬ 
apertures  (to  ensure  that  the  phase  variations  due  to  turbulence  over  each  sub-aperture 
are  small)  will  result  in  a  small  modelling  error,  at  low  light  levels,  this  will  result  in  a 
small  signal  in  each  sub-aperture  and  unacceptably  high  errors  in  the  slope  measurements. 
Alternatively,  reducing  the  number  of  sub-apertures  to  a  small  number  will  reduce  the 
noise  on  the  slope  measurements  but  give  a  larger  modelling  error.  For  a  given  signal-noise 
scenario,  there  will  therefore  be  an  optimum  choice  which  minimises  the  overall  effect  on 
the  average  mean-square  error  of  the  wavefront  estimate.  Attempts  to  study  this  problem 
have  been  made  but  the  analysis  is  somewhat  idealised  and  crucially  neglects  the  effect  of 
read-noise  for  CCD  devices  which  are  fast  becoming  the  de  facto  standard.  In  the  following 
article,  we  present  results  of  a  study  to  determine  optimal  Hartmann  sub-aperture  size  for 
different  signal-noise  levels.  This  is  based  on  the  use  of  a  high-fidelity  computer  model 
which  incorporates  the  effects  of  atmospheric  turbulence,  shot  noise  and  read-noise. 
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Abstract 

A  high-fidelity  computer  model  has  been  used  to  investigate  the  performance 
of  a  Shack-Hartmann  sensor  operating  under  low  light  level  conditions.  The  trade¬ 
off  between  sub-aperture  size  and  slope  noise  is  examined  in  detail  and  optimal 
sizes  are  determined  for  a  number  of  imaging  scenarios.  Our  results  are  discussed 
in  relation  to  related  literature. 
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1.  Introduction 


It  is  well  known  that  atmospheric  turbulence  severely  limits  the  performance  of 
ground  based  telescopes  and  there  now  exists  a  considerable  literature  on  both  tech¬ 
nological  and  computational  techniques  to  compensate  for  atmospheric  effects1 .  Ac¬ 
curate  estimation  of  atmospherically  perturbed  wavefronts  is  of  critical  importance 
(whether  directly  or  implicitly)  in  any  attempt  to  correct  incident  wavefronts  in  real 
time  using  adaptive  optics.  It  is  also  a  crucial  aspect  of  the  post-processing  tech¬ 
nique  of  deconvolution  from  wavefront  sensing  (also  termed  self-referenced  speckle 
holography)  originally  developed  in  reference2.  Although  wavefront  sensing  devices 
such  as  the  shearing  interferometer  and  increasingly  the  curvature  sensor3  have  been 
developed  and  used  for  astronomical  and  surveillance  imaging  systems,  the  practi¬ 
cal  and  conceptual  simplicity  of  the  Shack-Hartmann  device  and  it’s  competitive 
performance  still  make  it  a  first-choice  in  many  systems. 

In  this  paper,  we  use  a  high-fidelity  computer  model  employing  a  modal  wave- 
front  reconstructor  to  investigate  Hartmann  sensor  performance  at  low  signal  levels 
and  consider  the  corresponding  implications  for  optimal  design.  A  key  factor  in  em- 
ploving  Hartmann  sensors  for  adaptive  optics  applications  is  to  optimise  the  pupil- 
plane  sub-aperture  size.  When  the  signal  is  limited,  we  see  that  by  increasing  the 
Hartmann  sub-aperture  size,  the  signal-noise  of  the  individual  slope  measurements 
at  the  sub-apertures  is  increased  but  at  the  expense  of  larger  sampling  intervals  and 
hence  larger  modelling  error.  It  is  evident  that  a  balance  must  be  achieved  between 
these  two  competing  factors  if  the  best  performance  is  to  be  attained. 

Previous  authors  have  addressed  the  question  of  Hartmann  sub-aperture  size 
for  adaptive  optics  applications.  Gardner  et.  al.4  discuss  laser-guide  star  adaptive 
optics  systems  and  present  results  for  that  sub-aperture  size  which  will  minimise 
the  required  sub-aperture  photon  count  to  achieve  a  number  of  target  rms  wave- 
front  error  values.  Their  analysis,  however,  does  not  include  the  effects  of  read- 
noise.  Ellerbroek5  discusses  the  effect  of  varying  j-  (i.e.  varying  the  number  of 


sub-apertures  across  the  telescope  pupil  plane  but  this  effectively  only  studies  the 
fitting  error.  Gavel  et  al6  derive  analytic  results  for  optimum  sub-aperture  size  but 
under  the  restrictive  assumption  that  ?'o  <  1  and  that  CCD  read-noise  may  be  ne¬ 
glected.  They  also  present  results  for  the  optimum  sub-aperture  size  as  a  function 
of  laser  power  for  the  Keck  telescope  system. 

In  this  paper,  we  use  a  high-fidelity  computer  model  incorporating  turbulence, 
shot  and  CCD  read-noise  and  employing  a  modal  wavefront  reconstructor  to  inves¬ 
tigate  Hartmann  sensor  performance  at  low  signal  levels.  The  results  are  discussed 
with  regard  to  implications  for  optimal  design  of  systems  operating  in  low  signal- 
noise  regimes. 

2.  Computer  Model 

2.1  Generation  of  Hartmann  images 

Randomly  distorted  optical  wavefronts  obeying  the  Kolmogorov  power  spec¬ 
trum  were  generated'  and  considered  incident  at  the  pupil  plane  of  the  telescope. 
As  is  commonly  assumed  in  problems  of  imaging  through  turbulence,  amplitude 
fluctuations  were  neglected.  The  Hartmann  sub-aperture  images  were  calculated  by 
considering  the  appropriate  segment  of  a  high  (i.e.  infinite)  light  level  wavefront  and 
using  standard  Fourier  optics  to  generate  the  image  plane  intensities  as  - 

I((,v)ot\FT[A(x,y)}\2  (1) 

where  I(£,  77)  is  the  intensity  in  the  image  plane,  A(x,y)  is  the  complex  amplitude 
of  the  incident  wavefront  portion  and  FT  denotes  the  operation  of  Fourier  trans¬ 
formation. 

Photon  limited  images  are  generated  in  a  standard  way  from  the  infinite  fight 
level  images  via  generation  of  a  cumulative  distribution  function  which  is  then 
randomly  sampled  in  order  to  determine  the  location  of  each  photoelectron  event  in 
the  image  plane.  Poisson  sampling  for  the  number  of  detected  events  in  each  sub- 
aperture  is  also  included  in  the  model.  In  our  study,  we  have  assumed  a  CCD  camera 
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with  read-noise  of  1.5  electrons  r.m.s.  per  pixel  per  frame  -  this  figure  corresponds 
to  an  upper  limit,  on  the  performance  of  current  state-of-the-art  devices8.  A  number 
of  statistical  models  have  been  described  for  read  noise  in  CCD  arrays8-9.  In  our 
study,  read-noise  is  simulated  by  adding  an  appropriately  scaled  Gaussian  deviate 
to  each  pixel  in  the  CCD  image  frame. 

Calculation  of  the  sub-aperture  centroid  positions  (C(£,  r/)  —  (C^,Cn))  was 
performed  using  the  discretised  version  of  the  classical  centre  of  mass/first  moment. 

r  E€E,e*(^) 

EelVK,*) 

r  _  Eg  XV? jfog)  (0) 

""  EfE^K.v)  vw 

Initially  centroiding  was  performed  over  a  comparatively  large  area  to  locate 
the  approximate  position  of  the  centroid.  The  centroiding  envelope  was  then  re¬ 
duced  in  size  and  centroiding  was  repeated  centred  on  the  newly  estimated  spot 
location.  This  process  was  then  repeated  iteratively  until  the  en\ elope  had  reached 
the  optimal  dimension  1 ' '  and  the  centroid  location  had  stabilised.  This  typically 
takes  between  4  and  7  iterations  depending  on  the  stability  tolerances  specified  and 
the  rate  at  which  the  centroiding  envelope  is  reduced. 

Conversion  of  the  centroid  estimate  to  a  local,  average  wavefront  slope  is 
straightforward  and  the  slope  data  is  then  available  to  form  an  estimate  of  the 
incident  wavefront. 

2.2  Modal  wavefront  estimation 

A  substantial  body  of  work  has  now  appeared  in  the  literature  on  the  subject  of 
wavefront  estimation  and,  beginning  with  Wallner's  paper,  specifically  on  optimal 
wavefront  estimation5’13-16.  Within  this  body  of  literature,  due  to  variations  in 
notation  and  approach  (i.e.  integral  or  linear  algebraic  formulation),  to  the  wavefront 
model  (whether  zonal  or  modal)  and  the  specific  meaning  attached  to  optimal,  the 
common  thread  amongst  them  is  not  apparent.  For  clarity,  we  here  present  what  is, 
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to  the  best  of  our  knowledge,  the  most  direct  and  simple  derivation  of  a  maximum 
a  posteriori  (MAP)  estimator  for  a  linear  model.  Other  common  solutions  used  in 
wavefront  estimation  may  be  considered  either  as  alternative  forms  or  as  limiting 
cases  of  this  estimator. 

In  the  so-called  modal  approach,  the  random  phase  function  is  approximated 
by  an  expansion  of  the  form  - 

N 

*($)  =  (3) 

Jfe  =  l 

where  <&(x)  is  the  atmospherically  distorted  wavefront  and  {Pk{x)}  are  the  chosen 
basis  functions  and  the  system  is  modelled  as  - 

Ax  +  e  =  m  (4) 

where  x  is  the  vector  of  modal  coefficients,  A  is  the  system  or  design  matrix,  m  is 
the  vector  of  phase  gradients  and  e  is  the  error  vector  on  the  slopes. 

We  make  the  following  assumptions  on  this  linear  model  (all  of  which  are 
reasonable  for  Hartmann  and  shearing  sensors)  - 

•  Signal  and  noise  vectors  are  both  zero  mean  -  {x)  =  (e)  —  0. 

•  Noise  and  signal  are  uncorrelated  (xT e)  =  0. 

•  Noise  covariance  V  —  (eeT)  is  assumed  known. 

•  A  priori  covariance  of  the  unknown  vector  Cx  =  (xxT)  is  assumed  known.  This 
covariance  was  explicitly  calculated  for  the  Zernike  basis  and  Kolmogorov  spectrum 
by  Noll11.  A  suitable  unitary  transformation  enables  it  to  be  calculated  for  any 
other  orthogonal  basis  (e.g.  the  Karhunen-Loeve  basis)'. 

Accordingly,  we  seek  a  linear  estimator  x  =  Lm  which  will  minimise  the  error 
covariance  matrix  P  defined  by  - 

P  =  ((. x-x)(x-x)T)  (5) 

We  substitute  x  —  L(Ax  +  e)  in  eq.  5  to  obtain  - 

P  =[I  —  LA\CX[I  -  LA}t  +  LVLt 
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Taking  variations  in  matrix  P  with  respect  to  L  and  setting  to  zero14  .  w'e  have  - 
SP  -  SL(ACX[I  -  LAf  +  VLt)  +  ([/  -  LA\CxAt  +  LV)6LT  =  0  (7) 

which  yields  an  optimal  matrix  L  - 

Lmap  =  L  —  CxAt[ACxAt  +  V']-1  (7) 

Substitution  of  the  matrix  Lmap  into  eq.  5  defining  P  yields  an  error  covariance 
matrix  given  by  - 

Pmap  =  CX[I  -  At [AC xAt  +  V\~lAC,\  (8) 

It  is  easy  to  show  (see  appendix)  that  two  equivalent  forms  for  Lmap  and 

Pmap  are  - 

X  =  LMApm  =  [AtV~1A  +  Cp]-1ATV~1m  (9) 

P  =  Pmap  =  [ATV~1A  +  Cj1]"1  (9) 

Wallner’s  lengthy  analysis  reduces  to  this  estimator  under  the  same  assump¬ 
tions  (see  16)  and  this  may  also  be  derived  directly  from  Bayes  theorem14'15.  Fur¬ 
ther.  since  since  both  the  prior  and  likelihood  laws  on  the  modal  coefficients  may 
be  well-modelled  as  gaussians11,  the  choice  of  a  different  Bayesian  criterion  (such  as 
finding  the  mean  or  median  of  the  posterior  distribution  rather  than  the  mode)15 
is  academic  and  produces  the  same  result. 

In  the  limit  that  the  a  priori  covariance  tends  to  infinity,  the  MAP  estimator 
tends  to  the  best  linear  unbiassed  estimator  (BLUE)18  in  the  mean-squared  sense  - 

x  =  [ATV“lA]~lATV-1m  (10) 

If  we  assume  that  the  noise  on  each  measurement  is  uncorrelated  and  the  noise 
associated  with  each  measurement  is  equal,  we  have  V  =  a 2 1  and  eq.  10  gives  - 

x  =  [AT  A}-1  ATm  (11) 
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which  is  recognisable  as  the  standard,  least-squares  solution.  The  standard  least- 
squares  solution  is  thus  the  limiting  form  of  the  MAP  estimator  when  no  a  priori 
knowledge  of  the  parameter  statistics  and  noise  is  imposed  or  known. 

Both  the  BLUE  and  MAP  estimators  were  implemented  in  our  computer  model 
using  the  Zernike  and  Karhunen-Loeve  bases. 

3,  Simulation  Results 

In  practice,  Hartmann  wavefront  sensing  systems  are  limited  to  a  finite  number 
of  sub-apertures  and  a  finite  number  of  modes  only  can  be  estimated.  We  conducted 
a  number  of  studies  simulating  an  8x8  Hartmann  sensor  having  52  active  subaper¬ 
tures.  The  primary  goal  was  to  study  the  variation  in  maximum  possible  compen¬ 
sated  Strehl  ratio  as  a  function  of  flux  and  sub-aperture  size.  The  studies  were 
conducted  using  the  standard  least-squares  modal  reconstructor.  As  this  estimator 
does  not  include  a  priori  knowledge,  no  constraints  are  placed  on  the  modal  coef¬ 
ficients.  It  is  therefore  important  in  this  case  to  avoid  modelling  errors  which  may 
be  due  to  the  use  of  either  too  few  or  too  many  modes.  Accordingly,  we  first  study 
modelling  errors  for  our  system.  In  fig.  1  we  show  the  maximum  compensated  Strehl 
ratio  as  a  function  of  the  number  of  estimated  modes  with  sub-aperture  size  fixed 
=  r0.  The  Strehl  ratios  given  were  generated  from  the  residual  wavefront  (known 
wavefront  subtracted  from  the  estimated)  and  do  not  include  the  effects  of  anj'  real 
phase  compensating  system.  The  incident  wavefronts  correspond  to  high  light  level 
and  obev  the  Kolmogorov  power  spectrum  of  phase  fluctuations.  The  wavefront 
reconstruction  was  performed  using  the  standard  least-square  estimatoi  and  each 
data  point  is  averaged  over  400  realisations.  We  observe  a  peak  at  approximately 
45  modes  after  which  over-fitting  effects  increase  the  errors  significantly.  We  also 
note  that  although  the  peak  occurs  for  approximately  the  same  number  of  modes, 
there  is  a  significant  difference  between  the  Karhunen-Loeve  basis  and  the  Zeinike 
polynomials. 

Fig.  2  shows  the  effect  of  variations  in  the  photoelectron  flux  on  the  optimal 


value  of  l/r0.  Here,  an  8  x  8  Hartmann  sensor  having  52  active  subapertures  was 
used.  Read-noise  was  simulated  at  1.5  electrons  r.m.s.  per  pixel  per  frame.  As  intu¬ 
itively  expected,  the  peak  moves  towards  higher  values  of  l/r0  as  the  flux  falls  and 
lower  values  as  the  flux  increases.  It  is  clear  from  figure  2  that  the  correct  choice 
of  sub-aperture  size  l /vq  is  of  critical  importance  to  any  attempt  to  accurately  es¬ 
timate  and  compensate  wavefronts  under  conditions  of  low  incident  flux.  We  also 
note  that  the  pronounced  skew  in  the  figures  indicates  that  choosing  —  stna/Zerthan 
the  optimal  value  has  more  serious  consequences  than  slight  over-estimation  of  the 
optimum  value.  Fig. 3  shows  the  optimal  values  of  l f  vq  as  a  function  of  flux  levels 
for  a  number  of  CCD  read  noise  values.  The  read-noise  in  fig.  3  is  given  in  electrons 
per  pixel  per  frame.  In  any  particular  imaging  scenario,  the  temporal  bandwidth 
required  and  the  luminosity  of  the  guide  star  will  determine  the  frame  exposure 
time  and  hence  the  flux.  The  read-noise  of  a  CCD  device  for  the  chosen  integration 
time  can  then  be  determined.  In  this  way,  fig.  3  may  be  used  to  directly'  estimate 
or  extrapolate  so  as  determine  the  optimal  sub-aperture  size  for  the  given  set  of 
conditions.  We  point  out  that  the  Strehl  ratios  indicated  on  the  graphs  are  limit¬ 
ing  values  for  real  adaptive  optics  systems  which  will  not  be  achievable  in  practice 
since  thev  entirely  neglect  finite  bandwidth  effects  and  limitations  in  the  deformable 
mirror /correction  device.  However,  the  trend  observed  as  l/v o  increases  is  of  direct 
relevance  to  reed  AO  systems. 

Since  the  optimal  sub-aperture  size  clearly  depends  on  sub-aperture  signal- 
noise.  it  follows  that  the  optimal  values  indicated  in  figure  3  are  dependent  on  the 
specific  centroiding  algorithm  and  the  reconstructor  used.  We  chose  to  use  the  least- 
squares  wavefront  reconstructor  because,  despite  the  significant  theoretical  studies 
that  have  been  made  on  optimal  reconstructors,  implementation  of  realistic  and 
reliable  a  prion  knowledge  is  quite  complex  and  with  notable  exceptions’  has  rarely 
been  emploved  in  practical  sy'stems.  We  note  that  superior  slope  estimation  based  on 
matched  filtering  techniques  rather  than  centre  of  mass  should  yield  higher  signal- 
noise  at  low  light  levels.  However,  since  the  trade-off  between  sub-aperture  size  and 
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slope  estimation  noise  is  such  that  less  slope  measurement  noise  allows  smaller  sub¬ 
aperture  sizes,  we  predict  that  the  main  effect  of  such  methods  will  be  simply  to 
shift  the  optimal  sub-aperture  size  to  slightly  smaller  values  than  those  indicated. 
The  precise  amount  will  depend  on  the  relative  performance  of  the  matched  filter 
but  can  reasonably  be  expected  to  be  small. 

4.  Summary 

We  have  shown  that  the  optimal  sub-aperture  size  in  a  Hartmann  sensor  op¬ 
erating  at  low  fight  levels  is  critically  dependent  on  the  prevailing  signal  to  noise 
level  and  have  given  the  optimal  value  for  a  number  of  situations.  Critically,  our 
results  include  the  effect  of  CCD  read-noise.  The  optimum  value  for  a  number  of 
typical  imaging  scenarios  may  be  estimated  empirically  from  the  results  presented 
in  fig. 3-  This  may  be  >  3ro  in  very  low  signal-noise  conditions.  We  note  that  the 
specific  results  we  have  obtained  in  fig.3  for  the  optimal  sub-aperture  size  for  the 
zero  read-noise  case  are  in  good  agreement  with  those  obtained  by  Gardner  et  al 
(table  1  in  4). 
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FIGURE  CAPTIONS 


Figure  1 

Fitting  error:  The  maximum  possible  compensated  Strehl  ratio  with  number  of 
basis  functions.  The  sensor  has  52  active  sub-apertures  of  size  r0.  A  least-squares 
reconstructor  was  used. 

Figure  2 

Maximum  possible  compensated  Strehl  ratio  with  pupil-plane  subaperture  size 
for  signal  levels  of  10,  15,  20,  25  and  30  photoelectrons  per  rg.  Read-noise  of  2 
electrons  r.m.s  per  pixel  per  frame  was  generated. 

Figure  3 

Optimum  subaperture  size  as  a  function  of  photoelectron  flux  for  varying 
amounts  of  read-noise.  Flux  is  in  photoelectrons  per  and  read-noise  in  electrons 
r.m.s  per  pixel  per  frame. 
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Compensated  Strehl  Ratio 


2.  GRAM-CHARLIER  EXPANSION  FOR 
WAVEFRONT  SLOPE  ESTIMATION 

The  primary  estimation  task  in  Hartmann  sensing  is  that  of  finding  the  barycentre  (or 
"centroid" )  of  the  sub-aperture  point-spread  functions.  The  theory  of  Hartmann  sensing 
shows  that  the  vector  shift  in  the  position  of  this  barycentre  from  that  obtained  with  a 
known  reference  wavefront  is  linearly  proportional  to  the  average  wavefront  slope  over  the 
sub-aperture.  Errors  in  the  estimate  of  these  slopes  will  propagate  in  a  known  way  into 
the  estimate  of  the  wavefront  itself  and  it  is  thus  desirable  that  they  should  be  kept  as 
small  as  possible. 

The  conventional  method  for  calculating  the  barycentre  is  to  take  a  first  moment  or 
centre  of  gravity  calculation  on  the  data.  Speaking  generally,  one  of  the  advantages  of 
this  approach  is  that  it  makes  no  assumptions  on  the  nature  of  the  underlying  distribution 
from  which  the  data  has  been  generated.  However,  it  has  problems  associated  with  the 
introduction  of  a  systematic  bias  and  is  sensitive  to  noise.  Hartmann  sensors  operating 
with  pupil-plane  sub-apertures  sizes  of  a  few  r0  or  less  will  generally  produce  a  peaked 
PSF  whose  approximate  shape  may  be  predicted  quite  well.  There  is  thus  some  inherent 
a  priori  information  available  about  the  underlying  distribution.  When  signal  levels  are 
high,  there  is  little  advantage  to  be  gained  by  seeking  an  alternative  to  the  first  moment 
estimator.  However,  in  low  signal-noise  regimes,  the  sensitivity  to  noise  becomes  manifest 
and  an  estimator  with  better  noise  suppression  is  desirable.  The  following  article  presents 
results  obtained  with  a  matched  filter  based  on  modelling  the  PSF  as  a  Gram-Charlier 
expansion.  The  performance  is  compared  with  the  first  moment  estimator  for  a  number  of 
signal- noise  scenarios. 
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A  study  has  been  made  of  a  Gram-Charlier  matched  filter  for  Shack-Hartman 
sensing  of  wavefront  slopes.  The  method  is  based  on  modelling  the  Point  Spread 
Function  (PSF)  by  an  expansion  in  terms  of  Gauss-Hermite  polynomials.  We 
present  results  for  several  sub-aperture  /  coherence  area  sizes  both  with  and  with¬ 
out  CCD  read  noise.  More  accurate  estimation  of  the  local  slopes  may  be  achieved 
at  low  light  level  as  compared  to  the  standard  first  moment  estimator. 
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sion. 

High  quality  performance  in  Hartmann  sensors  is  dependent  on  the  accu¬ 
rate  estimation  of  wavefront  slopes  and  thus  on  the  sub-aperture  centroids. 
This  is  often  calculated  satisfactorily  by  performing  a  first  moment  calcula¬ 
tion  on  the  sub-aperture  intensity  distribution.  However,  at  low  flux  levels 
the  first  moment  estimator  for  the  centroid  has  poor  noise  propagation  char¬ 
acteristics  [1].  Essentially  this  originates  from  the  fact  that  no  assumptions 
or  constraints  are  placed  on  the  underlying  distribution  from  which  the  data 
has  been  sampled.  In  noisy  scenarios,  use  of  the  first  moment  estimator  will 
mean  that  pixels  far  from  the  underlying  peak  will  make  a  large  contribution 
to  the  moment  but,  as  the  underlying  distribution  is  generally  unimodal, 
these  pixels  will  have  low  signal-noise  ratios.  This  approach  to  estimation 
also  has  problems  associated  with  the  introduction  of  a  systematic  bias  [2,3]. 
There  is  thus  a  clear  motivation  for  seeking  a  matched  filter  which  effectively 
weights  the  data  so  as  to  achieve  better  noise  propagation  characteristics  in 
noisy  conditions.  This  fact  has  been  recognised  by  workers  in  a  variety  of 
fields  including  nuclear  experimentation,  position  tracking,  microscopy  and 
fluid  mechanics  [4-7]. 

For  adaptive  optics  applications  using  Hartmann  sensors,  approaches  using 
artificial  neural  networks  [8]  and  Maximum  a  Posteriori  (MAP)  estimation 
of  Hartmann  centroids  [9]  have  recently  been  proposed.  In  this  paper  we 
consider  a  Shack-Hartmann  sensor  utilising  a  CCD  array  as  a  detector  and 
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investigate  an  alternative  approach  to  estimating  the  subaperture  centroid. 
This  is  based  on  the  modelling  of  the  sub-aperture  PSF  as  a  truncated  Gram- 
Charlier  expansion  [10]. 

We  consider  the  light  focussed  onto  the  CCD  detector  from  one  sub¬ 
aperture  of  the  Shack-Hartmann  sensor.  We  will  assume  that  this  photon 
distribution  be  represented  over  a  square  matrix  of  N 2  pixels,  with  each  el¬ 
ement  registering  a  number  of  detected  photoelectrons.  CCD  read  noise  is 
simulated  by  adding  a  Poisson  sampled  random  deviate  to  each  element  of 
the  matrix.  The  spatial  distribution  of  the  photoelectrons  I(xi,yj)  is  then 
given  by  - 


I(xi,yj)  =  h(xi,yj)  +  N(xi,yj)  (i  =  1, ..,  N;j  —  1,  ..N)  (1) 

Where  N(xi,yj)  is  the  Poisson  random  value  having  probability  density  - 

P(k)  = 

The  parameter  m  specifies  the  mean  (and  the  standard  deviation)  of  the 
noise  process. 

The  centroid  of  the  distribution  is  given  by  the  first  moment  estimator: 


cx 


ZhxJixi.yj)  Ef=i  vJ{xi,  yj) 
ZiLi  l{xu yj)  'Cy  Thlixi.Vj)  ' 


(2) 
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The  Gram-Charlier  expansion  is  a  convenient  method  for  representing  near 
gaussian  distributions  [10].  We  note  that  when  Hartmann  wavefront  sensors 
are  used  in  problems  of  adaptive  optics  and  imaging  through  atmospheric 
turbulence,  it  is  common  to  employ  them  such  that  the  sub-aperture  size 
(scaled  to  the  telescope  pupil  plane)  will  generally  range  from  approximately 
1  to  3 r0  where  r0  is  the  Fried  coherence  parameter.  For  sub-aperture  sizes 
in  this  range,  the  Hartmann  PSFs  will  generally  be  peaked  and  unimodal 
and  may  reasonably  be  considered  as  approximately  gaussian  in  form.  In 
our  study,  we  have  used  a  two-dimensional  Gram-Charlier  expansion  of  the 
intensity  distribution  I(x,  y )  in  the  CCD  image  plane  of  the  following  form  - 


I(x.  y )  =  Go(x,y)  <  1  +  ^ 


,afkAj’k(x,y) 


— TT - \\  r  (3) 

Where  j  and  k  are  all  possible  coombinations  such  that:  j  +  k  —  i  and 
0  <  j  <  i,  0  <  k  <  i.  and  where  G0  represents  the  best  least-squares 

gaussian  fit  to  the  distribution  I(x.y): 


G0(x.y)  =  ae 


(j) 


The  free  parameters  a,  u.v  ox  and  ay  for  the  gaussian  fit  are  found  by 
minimising  a  standard  y2  error  metric  - 

X2{<y,  U.  V,  <jx,  oy)  =  X)  [I(xii  Vi)  -  G0{xi,  yj)]2  (5) 

2=1  j—  1 

The  A?k  functions  are  defined  as  - 


Afk (x,  y)  =  Hj(x)Hk{y):  j  +  k  =  i,  0  <  j  <  i,  0  <  A:  <  i  (6) 
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where  the  H i  represent  the  Hermite  polynomials.  The  Hermite  polyno¬ 
mials  are  a  complete  orthogonal  set  over  the  real  axis  with  respect  to  the 
weight  function  /\/27r.  Thus. 

<W,Hk>  =  4=  /  e-^Hj[x)Hk(x)dx 

v27T  JR 

=  j'-tijk  (7) 


The  Aj’k  are  two  dimensional  polynomials  of  degree  i— j+k  and  are  a 
complete  orthogonal  set  over  the  space  of  all  real  number  pairs,  R[X,Y]  - 


<  Aik,A{,’k' >  =  /R  JR  ^~-HJ{.r)Hf{x)Hk{y)Hkl{y)dxdy 


Ir  2tt 

=  j\k\8(j-j')S(k-k') 


(8) 


The  coefficients  in  the  expansion  described  by  eq.3  are  thus  given  by 
evaluating  the  integral  - 


/  I(x,  y) 

\Go(.'C,  y) 


1 )  A\k(x.  y)dxdy 


which  in  discrete  form  becomes  - 


(9) 


N  N 


VGo  (xi,yP) 


—  ( - 1)  Mk(xi,yP) 


(10) 


A  computer  simulation  was  performed  to  compare  the  Gram-Charlier  matched 


filter  and  first  moment  methods  for  estimating  the  sub-aperture  PSF  cen¬ 
troid.  Wavefronts  were  generated  having  statistics  that  obey  the  Kolmogorov 
model  [11]  and  the  resulting  reference  (i.e.  infinite  light  level)  sub-aperture 
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PSFs  were  generated.  The  precise  centroid  of  these  reference  PSFs  was  then 
calculated.  Taking  the  infinite  light  level  PSF  as  an  unnormalised  probability 
density,  the  PSFs  were  then  photon-sampled  using  standard  techniques  and 
CCD  read  noise  was  added  to  produce  a  sub-aperture  photon  distribution. 
The  centroid  of  this  distribution  was  then  calculated  using  both  the  first 
moment  estimator  and  the  Gram-Charlier  method  and  compared  against  the 
exact  value  calculated  from  the  reference  PSF.  A  comparison  was  made  for 
three  sub-aperture  size  ratios  of  ^  =  1,2  and  3.  A  mean  signal  level  of  30 
photoelectrons  per  sub-aperture  was  used  in  all  cases.  Read  noise  was  sim¬ 
ulated  using  a  Poisson  model  at  a  level  of  1.5  electrons  r.m.s  per  pixel  per 
frame.  The  simulation  was  configured  such  that  the  generated  PSFs  were 
contained  within  an  area  of  3*3  CCD  pixels  [12].  This  area  was  located  and 
the  first  moment  estimation  of  the  centroid  was  performed. 


•  Figure  1  compares  the  two  estimators  for  the  situation  in  which  no  detec¬ 
tor  read  noise  is  added  to  the  sampled  data.  The  data  points  corresponding 
to  the  ordinate  labelled  -1  represent  the  rms  error  in  CCD  pixels  of  the 
centroid  from  use  of  the  first  moment  estimator  and  the  data  points  at  the 
ordinate  0  represent  the  rms  error  resulting  from  a  least  squares  Gaussian  fit 
to  the  data.  Data  points  at  subsequent  integer  values  of  the  ordinate  give 
the  rms  error  as  the  radial  order  of  the  Gauss-Hermite  fit  is  increased.  The 
graph  shows  an  improvement  (~15%)  for  the  Gaussian  and  first  order  Her- 
mite  functions.  We  see  that  adding  more  terms  does  not  significantly  affect 
the  accuracy. 
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•  Figure  2  was  generated  under  similar  conditions  to  figure  1  but  compares 
the  two  estimators  when  simulated  CCD  read  noise  has  been  added  to  the 
photon  limited  data.  This  graph  confirms  that  the  gaussian  matched  filter 
yields  a  significant  improvement  in  the  slope  estimate.  Addition  of  the  first 
order  Hermite  terms  yields  a  further  minor  gain  after  which  the  higher  order 
terms  reduce  the  accuracy.  This  behaviour  is  intuitively  expected  since  the 
fitting  of  higher  order  terms  will  increasingly  tend  to  reproduce  the  noise.  It 
is  also  consistent  with  the  fact  that  their  contribution  was  negligible  in  figure 
1  where  no  read-noise  was  present. 

•  Figure  3  illustrates  the  estimation  procedure  graphically  using  contour 
plots.  A  high  light  level  Hartmann  subaperture  image  (first  row:  left)  is 
sampled  for  30  photoelectrons  and  CCD  read  noise  added  using  a  mean  level 
of  1.5  electrons  r.m.s  per  pixel  per  frame  (first  row:  middle).  The  result 
of  performing  a  gaussian  fit  (first  row:  right)  is  followed  bv  fits  including 
successively  higher  order  Hermite  terms. 

•  Table  1  compares  some  corresponding  r.m.s  wavefront  errors  for  the  first 
moment  estimator  and  for  the  Gram-Charlier  method.  Results  are  given  for 
the  least-square  and  MAP  reconstructors  [14.15]  for  a  number  of  different 
operating  conditions.  We  see  that  as  the  signal  decreases,  the  gain  of  the 
Gram-Charlier  estimator  over  the  first  moment  estimator  increases. 

We  have  investigated  a  matched  filter  for  Shack-Hartinann  slope  estima¬ 
tion  at  low  light  levels  based  on  a  truncated  Gram-Charlier  expansion.  The 
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orthogonality  of  the  functions  allows  independent  and  direct  evaluation  of  the 
coefficients  by  integration.  Our  results  indicate  that  a  significant  increase  in 
accuracy  over  the  standard  first  moment/centre  of  gravity  estimator  can  be 
made  for  subaperture  sizes  from  1  -  3r0.  A  recent  study  of  low-light  level 
wavefront  sensing  suggests  that  the  optimum  subaperture  size  will  generally 
fall  within  this  region  [13]. 

In  the  calculation  of  the  Gram-Charlier  coefficients  we  found  it  necessary 
to  restrict  the  region  of  integration  indicated  by  eq.  11  to  a  region  within 
approximated  2.5  standard  deviations  about  the  centre  of  the  Gaussian  fit. 
The  presence  of  noise  means  that  integration  over  a  larger  interval  results 
in  instability  in  the  Gram-Charlier  coefficients.  In  the  current  work,  the  ap¬ 
propriate  domain  of  integration  has  been  established  on  an  empirical  basis. 
We  note  that  Bayesian  MAP  estimation  techniques  have  been  successfully 
applied  both  to  centroid  and  to  wavefront  estimation  [9.14]  and  these  should 
provide  a  natural  way  of  allowing  the  domain  of  the  fitting  procedure  to  be 
extended  whilst  constraining  the  estimated  coefficients.  Such  an  approach 
requires  that  a  large  ensemble  of  high  light  level  subaperture  PSFs  be  mea¬ 
sured  to  generate  prior  knowledge  in  the  form  of  an  ensemble  covariance  of 
the  Gram-Charlier  coefficients. 
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Figurel:  Mean-square  error  in  centroid  estimate  using  the  first  moment 
estimator  and  for  different  orders  (0  to  6)  of  the  Gram-Charlier  fit.  The  flux 
level  is  30  photoelectrons  per  sub-aperture.  No  read-noise  is  present  and  ^ 
varies  from  1  to  3. 

Figure2:  Mean-square  error  in  centroid  estimate  using  the  first  moment 

estimator  and  for  different  orders  (0  to  6)  of  the  Gram-Charlier  fit.  The  flux 

level  is  30  photoelectrons  per  sub-aperture.  Poisson  read-noise  is  added  and 

—  varies  from  1  to  3. 
ro 

Figure3:  Graphical  illustration  of  the  effect  of  adding  more  Gram-Charlier 
terms:  Initially,  the  centroid  estimate  improves  but  the  addition  of  too  many 
terms  tends  to  reproduce  the  noise. 

Table  1:  R.m.s.  wavefront  reconstruction  error  for  the  first  moment  and 
Gram-Charlier  estimators.  Thee  results  are  calculated  for  two  wavefront  re¬ 
constructors:  least-squares  (zrms)  and  MAP  (maprms). 
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MS  ERROR  IN  CENTROID  POSITION  (RAD  /M^) 


Number  of  Flux  Corresponding  Visual  A  priori  phase  variance  R.m.s  Error  (rad)  R.m.s  Error  (rad) 

ubapertures  (p.e./i'o2  /frame)  Magnitude  (1.015  (D/r0)5/6  rad)  1st  moment  estimator  Gram-Charlier  estimator 

Exposure  time  =5  ms 


3.  BAYESIAN  ESTIMATION,  MODAL  PROJECTION  AND  SENSOR  GEOMETRY 

A  maximum  a  posteriori  (MAP)  Bayesian  wavefront  estimator  was  derived  in  section 
1  of  the  interim  report1.  The  error  covariance  matrix  of  this  estimator  was  derived  and 
evaluated  for  a  number  of  typical  Hartmann  imaging  scenarios  and  compared  to  the  best 
linear  unbiassed  estimator  (BLUE).  The  derivation  and  explicit  method  of  evaluation  is  not 
therefore  repeated  here.  Results  obtained  from  this  are  given  in  figs  3.1  -  3.3.  Significant 
gains  are  apparent  at  low  flux  levels.  These  results  do  not  include  the  effects  of  detector 
noise  (for  example  CCD  read  noise)  as  it  is  difficult  to  analytically  incorporate  this  into 
the  noise  model.  However,  the  basic  trend  of  these  results  which  indicates  a  significant 
improvement  at  low  fight  levels  can  be  expected  to  appear  when  detector  noise  is  present. 

The  superior  performance  of  the  MAP  estimator  is  dependent  on  having  the  correct 
form  for  the  a  priori  atmospheric  phase  covariance  matrix,  C.  Non- Kolmogorov  behaviour 
of  the  atmosphere  would  therefore  require  an  experimentally  determ  ined  form  for  the 
covariance  rather  than  a  theoretical/ analytical  one. 

The  issue  of  sensor  geometry  received  some  discussion  in  the  early  literature  [2-6]  where 
the  conclusion  was  reached  that  the  wavefront  error  is  fairly  insensitive  to  sensor  geometry 
(in  the  sense  of  size  and  placement  of  subapertures).  For  this  reason,  most  Haitmann 
sensors  use  the  simple  and  practical  “square-grid"  geometry  giving  a  uniform  sampling 
over  the  pupil.  Recent  work  [7,8]  has  suggested  however,  that  gains  can  potentially  be 
made  by  combining  a  Hartmann  geometry  conforming  to  the  Albrecht  cubature  nodes 
with  a  suitably  constructed  set  of  vector  polynomials  to  estimate  wavefront  modes  by 
direct  integration.  To  see  how  this  is  achieved  in  principle,  it  is  necessaiy  to  first  introduce 
a  set  of  vector  polynomials  which  can  restore  the  orthogonality  of  the  linear  problem  that 
arises  in  Hartmann  sensing. 

3.1  Vector  polynomials  orthogonal  to  the  gradient  of  the  basis  functions 

In  problems  of  wavefront  estimation,  it  is  common  to  expand  the  wavefront  in  terms 
of  a  set  of  orthogonal  functions.  The  most  common  of  these,  of  course,  are  the  Zernike 
circular  polynomials.  Speaking  generally,  orthogonal  bases  are  often  desirable  because 
they  allow-  the  expansion  coefficients  to  be  evaluated  by  simple  integration  of  a  product  of 
tw-o  functions  over  the  domain.  However,  since  Hartmann  sensors  (and  shearing  interfer¬ 
ometers)  provide  estimates  of  the  gradient  of  the  phase  rather  than  the  phase  itself,  the 
orthogonality  of,  for  example  the  Zernike  or  Karhunen-Loeve  basis  cannot  be  exploited. 
Modal  cross-coupling  can  occur  and  coefficients  must  be  obtained  by  solving  a  least-squares 
procedure. 

Consider  the  wavefront  to  be  expanded  in  terms  of  an  orthogonal  set  of  modes  - 

N 

$(:£)  =  ^2  OkPk(L) 

k  =  l 


(3.1) 


(3.2) 


Slope-based  wavefront  sensing  implies  that  our  model  is  thus 

N 


fc=i 


The  evaluation  of  modes  by  direct  integration  could  be  restored  if  we  can  derive  a  set 
of  auxiliary  functions  Ffx)  which  are  orthogonal  to  the  gradients  of  the  basis  functions. 
Let  these  functions  obey  the  relation  - 


V  •  Ffx)  =  Pfx) 


/here  the  orthogonality  relation  for  the  basis  is  - 

/  Pfx)Pj(x)d2x  =  Sij 

Jd 


(3-3) 


(3.4) 


and  D  denotes  the  domain  of  integration. 

Using  the  orthogonality  of  the  basis  in  eq.  3.1  and  substituting  eq.  3.3,  the  modal 
coefficients  are  given  by  - 


ai  =  [  <Kx)Pi(x)cPx  =  /  <?(z)V  •  Ffx)d2 
Jd  jd 


(3.5) 


This  may  be  written  as  - 

ai  =  J  V  •  (Ffx)  ■  <f>(x))d2x  -  V<f>{x)  ■  Fi(x)d2x 

and  applying  the  divergence  theorem,  we  obtain  - 

ai  =  j  <t>U)Fi{x)  -dl-  J  V<f>{x)  ■  Fi{x)d2x 

If  the  required  set  of  functions  Ff  x)  obeys 

V  •  Ff  x)  =  Pf  x) 

Fi(x)\n(C)  =  0  (3.6) 

where  Fi(x_)\n(C)  denotes  the  normal  component  of  the  function  to  the  closed  contoui  C 
of  the  integration  domain  D,  we  then  have  - 


a,  =  -  /  Vo(l 

JD 


)  ■  Fi(x)drx 


(3.7) 


and  we  may  evaluate  the  modes  by  direct  integration  as  required. 

There  are  in  fact  a  number  of  possible  sets  of  vector  polynomials  which  are  orthogonal 
to  the  gradient  of  the  Zernike  basis  [7,9,10].  This  is  because  there  is  no  unique  solution  to 
eq.  3.3  with  the  given  boundary  conditions  eq.  3.6.  Here,  we  make  use  of  the  set  originally 
derived  by  Gavrielides  and  refer  the  reader  to  the  original  source  for  the  derivation [9].  In 
the  appendix  A3,  we  prove  that  these  functions  are  irrotational  and  minimum  norm.  This 
property  of  the  functions  appears  to  have  been  unknown  until  now  but  is  important  because 
it  guarantees  that  of  all  possible  sets  of  vector  polynomials  orthogonal  to  the  gradient 
of  the  basis  functions,  they  will  have  the  smallest  noise  propagation.  For  completeness, 
the  explicit  form  for  generation  of  the  Gavrielides  functions  in  polar  coordinates  and  the 
corresponding  error  covariance  (noise  propagation)  for  the  modal  estimator  described  by 
eq.  3.7  are  also  given  in  the  appendices  A1  and  A2. 

3.2  Noise  propagation  and  Sensor  geometry 

By  making  use  of  the  Gavrielides  polynomials,  the  modal  coefficients  can  be  estimated 
by  direct  integration  over  the  pupil  according  to  eq.  3.7.  We  point  out  that  restoration  of 
the  orthogonality  via  Gavrielides  functions  allowrs  us  to  estimate  the  puie  modal  content 
of  a  distorted  wavefront.  In  other  words,  the  contribution  of  each  mode  is  estimated  inde¬ 
pendently  of  the  others.  For  a  non-orthogonal  problem,  least-squares  solutions  effectively 
allow  a  shuffling  of  modes  (in  a  manner  analogous  to  aberration  compensation  in  optical 
design)  so  as  to  achieve  the  minimum  norm  error  with  respect  to  the  data.  However,  as  far 
as  minimising  the  mean-square  wavefront  error  there  is  no  advantage  per  se  in  transform¬ 
ing  to  an  orthogonal  basis.  This  has  been  noted  by  Southwell  [9].  Any  advantage  must 
come  as  the  result  of  sampling  the  pupil  plane  in  a  way  which  gives  reduced  noise  on  the 
calculations  which  are  required  to  estimate  the  modes. 

In  Hartmann  sensing,  we  do  not  have  the  continuous  function  V$>(z)  but  rather  have 
it’s  value  at  a  number  of  sample  points.  The  simplest  approach  to  estimating  the  modes  us¬ 
ing  the  auxiliary  vector  functions  (and  a  reasonable  one  for  a  regular  Hartmann  geometry) 
is  to  approximate  the  integral  by  a  simple  two-dimensional  summation  - 

at  =  -  Y,  VttZj  )  •  Fi (xj  )dA  (3.8) 

j 

where  the  area  element  dA  is  simply  equal  to  the  area  of  the  pupil  (n  divided  by  the  total 
number  of  sample  points  (sub- apertures)  ). 

There  is  a  large  literature  devoted  to  efficient  and  accurate  numerical  integration  of 
functions  sampled  over  closed  domains.  In  general,  these  schemes  will  attempt  to  identify 
those  points  at  which  smooth  functions  (typically  polynomials)  need  to  be  sampled  within 
the  given  domain  in  order  that  they  may  be  integrated  most  accurately.  The  cubatures 


which  we  have  used  in  this  study  are  due  to  Albrecht  and  Engels  [11]  and  define  a  set  of 
x.y  points  lying  within  a  circular  domain  with  each  of  which  we  associate  a  weight  so  that 
we  approximate  the  integral  by  finite  sums  of  the  form 

N 

ak  =  -c^1  ^  wafaV<t){a)  •  Fk{a )  (3-9) 

<*=  1 

a  =  i  _>  jV  labels  the  points  at  which  the  integrand  is  evaluated  (known  as  the  nodes 
of  the  cubature),  wa  are  the  weights  associated  with  these  points  and  fa  is  1  if  cartesian 

coordinates  are  used  and  ra  if  polar  coordinates  are  used. 

In  the  appendix  A2,  the  mean-square  error  resulting  from  use  of  the  estimator  in  eq. 

3.8  is  derived.  It  is  given  in  matrix  form  by  - 

a^TrilF^F  +  C;1)-1}  (3-10) 

The  mean  square  error  resulting  from  use  of  the  Albrecht’s  grid  geometry  is  given  by 

a2  =  -cf  w2jlFk(a)  ■  Fk(a)  (3.11) 

O'—  1 

This  error  covariance  matrix  was  generated  using  Gavrielides  functions  and  both  the 
regular  square  Hartmann  geometries  and  Albrecht  s  grid  geometries  for  several  signal- 
noise/turbulence  scenarios.  Also  compared  were  the  errors  on  the  individual  modes  for  the 
MAP  bayesian  estimator  and  the  BLUE  estimator.  Our  studies  have  shown  that  although 
the  use  of  the  MAP  estimator  results  in  reduced  noise  propagation  (see  figs.  3.4  -  3.6)  for 
32,  43  and  80  subapertures  respectively,  the  geometrical  arrangement  of  subapertures  does 
not  appear  to  have  a  significant  effect  on  the  noise  propagation. 
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FIGURE  3.2 

r.m.s.  wavefront  error  against  photoelectron  signal  in  each  sub-aperture. 
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FIGURE  3.5 

Noise  propagation  on  the  individual  Zernike  modes  for  the  Gavrielides  functions 
with  regular  Hartmann  geometry. 
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FIGURE  3.6 

Noise  propagation  on  the  individual  Zernike  modes  for  the  GavrieUdes  functions 
with  regular  Hartmann  geometry.  The  Bayesian  MAP  solution  has  better  noise 
propagation  than  the  BLUE. 
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4:  SUMMARY 


4.1  Optimal  subaperture  size  for  Hartmann  sensors 

A  study  has  been  made  of  optimising  the  Hartmann  sub-aperture  size  for  wavefront 
sensing  applications  at  low  light  levels  and  submitted  for  publication  to  Optics  Communi¬ 
cations.  This  is  based  on  a  detailed  computer  model  simulating  the  effects  of  turbulence, 
Poisson  and  CCD  read  noise  on  the  wavefront  sensing  process.  Our  results  demonstrate 
that  optimal  wavefront  estimation  and  correction  is  critically  dependent  on  matching  the 
sub-aperture  size  to  the  prevailing  signal-noise  ratio  and  we  give  the  optimal  values  for  a 
number  of  typical  low  light  level  imaging  scenarios. 

4.2  Gram-Charlier  centroid  estimation 

A  study  has  been  made  of  a  Gram-Charlier  matched  filter  for  the  estimation  of  Shack- 
Hartmann  sub-aperture  barycentres  at  low  light  levels.  The  results  of  this  work  have 
been  submitted  and  accepted  for  publication  in  Optics  Letters.  A  significant  increase  in 
accuracy  over  the  standard  first  moment/centre  of  gravity  estimator  can  be  made  when 
the  subaperture  size  is  2  -  3r0.  The  residual  mean-square  error  in  the  wavefront  estimate 
has  been  calculated  and  the  limiting  visual  magnitude  at  which  this  performance  may  be 
achieved  has  been  calculated  for  a  number  of  imaging  conditions. 

4.3  Bayesian  and  modal  projection  estimation 

We  have  derived  and  implemented  a  Bayesian  estimator  for  reconstruction  of  atmo- 
sphericallv  perturbed  wavefronts.  The  prior  knowledge  is  implemented  in  this  estimator  by 
assuming  that  the  atmosphere  induces  phase  fluctuations  obeying  the  Kolmogorov  power 
spectrum.  The  covariance  matrix  for  the  atmospheric  modes  thus  obeys  the  form  originally 
derived  by  Noll  [1].  Calculations  of  the  error  covariance  matrix  for  a  number  of  different 
sub- aperture  sizes  £  and  a  number  of  different  turbulence  strengths  have  indicated  that 
the  Bayesian  estimator  has  a  reduced  mean-square  error  as  compared  to  the  standard 
least-squares  estimator.  This  improvement  is  particularly  significant  at  low  signal  levels. 
These  results  are  presented  and  discussed  fully  in  the  interim  report1. 

We  have  also  investigated  a  method  of  modal  projection  estimation  based  on  direct 
integration  using  a  set  of  auxiliary  vector  polynomials  which  are  orthonormal  to  the  gra¬ 
dients  of  the  Zernike  polynomials.  This  is  a  quick  and  convenient  procedure  which  ensures 
that  the  modes  are  estimated  individually  without  modal  cross-coupling.  However,  since 
adaptive  optics  applications  are  generally  concerned  with  minimising  the  global  mean 
square  wavefront  error  (rather  than  errors  on  particular  individual  modes)  we  have  ex¬ 
amined  the  overall  noise  propagation  of  this  method  when  combined  with  the  Albrecht’s 


cubatures  defined  on  the  unit  radius  circle.  We  have  developed  a  Bayesian  estimator  using 
the  auxiliary  vector  polynomials  and  have  shown  that  the  noise  propagation  is  superior  to 
that  obtained  by  direct  modal  projection  in  which  no  a  priori  knowledge  is  incorporated. 
Although  there  is  some  computational  advantage  in  using  the  auxiliary  vector  polynomi¬ 
als,  our  study  suggests  that  the  use  of  tailored  Hartmann  sub-aperture  geometries  has  no 
significant  effect  on  noise  propagation  and  will  not  lead  to  reduced  mean-square  error  in 
wavefront  estimation  as  compared  to  the  regular  Hartmann  geometry. 


APPENDICES 

Al  :  GAVRIELIDES’  POLYNOMIALS 

These  vector  polynomials  are  orthonormal  with  respect  to  the  gradient 
of  the  Zernike  polynomials.  They  are  defined  in  polar  coordinates  as  - 


1.  m/0 


1  t - -  _  x  I  cosmO  ifieven 

Fp  =  —  \/2(n  +  1)T„  (r)  <  s^nmQ  i fiodd 


1  / - -  ,  f  —msinmO  ifieven 

F'e  —  ~V2(n  +  (T )  |  rncosmO  ifiodd 


where: 


1  Tlf  (7.™(sl  In  -  2s  +  2)  frm~~1  -  r"-2s+1] 
Tn  ^  ~~  4  5  -  s  +  l)  -s  +  l) 


1  c™($)  [(n.  -  2s  +  2 )rm~1  -  mrn  2s+1] 

^  W  =  IS  5  («±i-»  +  l)  + 


2.  m=0 


.  1  , - -X  (-1  )5(n-s)!r"-25+1 

F  =  -  V«  +  1  zZ  — ^ - TT~  77 

p  *■  (f-s)l  (n  —  2s  +  2) 


Fi  =  o 


A2:  NOISE  PROPAGATION  FOR  GAVRIELIDES’  FUNC¬ 
TIONS 

Our  model  for  the  wavefront  is  a  finite  expansion  in  terms  of  a  known  set 
of  basis  functions.  Accordingly,  we  begin  from  the  discretised  linear  system 
given  in  eq.  1.4  - 

Ar  +  e  =  m  (5) 

where  x  is  the  hi  by  1  vector  of  unknown  modal  coefficients,  A  is  the  system 
or  design  matrix  of  size  M  by  N  and  containing  the  gradients  of  the  basis 
functions  evaluated  at  the  sensor  positions,  m  is  the  M  by  1  phase  gradient 
measurement  vector  and  e  is  the  error  vector.  We  make  identical  assump¬ 
tions  on  this  model  to  those  outlined  in  section  1.  Let  us  apply  a  linear 

transformation  to  this  system  of  equations  by  applying  a  matrix  FT  from 
the  left  on  both  sides  - 


FtAx  +  FTe  =  FTm 

The  result  is  also  a  matrix-vector  equation  which  we  denote  - 


Ax  +  e  =  m 


(6) 

(7) 


So,  we  seek  a  linear  estimator  x  =  Lm  which  will  minimise  the  error 
covariance  matrix  P  defined  by  - 

P  =  ((x-x)(: x-xf)  (8) 

We  substitute  x  =  L(  Ax  +  e)  in  eq.  8  to  obtain  - 

P  =  [/  -  LA]CX[I  -  LA?  +  LVLt  (9) 

where  V  —  (eeT)  =  FTVF  Taking  variations  in  matrix  P  with  respect  to  L 
and  setting  to  zero,  we  have  - 

6P  =  SL(ACX{I  -  LA?  +  VLT)  +  ([/  -  LA]CXAT  +  LV)6LT  =  0  (10) 

which  yields  an  optimal  matrix  L  - 

Lmap  ~  L  =  CxAt[ACxAt  +  V]  1 


(11) 


Substitution  of  x  =  Lmap  m  =  Lmap(Ax  +  e)  into  eq.4  defining  P  yields 
an  error  covariance  matrix  given  by  - 

Pmap  =  cx[l  -  At[ACxAt  +  VY'ACr)  (12) 

In  appendix  A,  we  show  the  equivalence  of  eqs  11  and  12  to  the  following 
two  forms  - 

x  =  Lmap  m  =  [AJ'V-1  A  +  Cx  *]  1ATV  1rh  (13) 

P  =  Pmap  =  [AtV~'A  +  C;1]-1  (14) 

Now  suppose  that  we  are  able  to  construct  a  set  of  functions  (represented 

in  sampled  form  by  the  columns  of  the  matrix  F)  which  are  orthogonal  to 
the  gradient  of  the  basis  functions  (i.e.  the  Gavrielides  functions).  In  other 
words,  we  are  able  to  find  a  matrix  such  that  FT A  =  L  the  identity  matrix. 
Since  we  derived  equations  13  and  14  for  a  quite  general  matrix  F,  we  may 
then  simply  substitute  back  A  =  FT  A,  e  =  FTt  and  m  =  FTm  into  eqs  13 
and  14  to  obtain  an  estimator  for  the  modal  coefficients  - 

x  —  FTm  (15) 

with  an  associated  error  covariance  of  - 

P  =  Pmap  =  [{FTVF)~1  +  C;1]'1  (16) 

The  noise  propagation  on  the  kth  mode  will  thus  be  given  by  the  A th  diagonal 
element  of  P  whilst  the  overall  mean  square  error  (excluding  under-modelling 
error)  will  be  given  by  the  trace  of  matrix  P.  ^ 

If  we  assume  that  the  errors  are  uncorrelated  and  equal  such  that  \  —  a  I 
and  the  a  priori  covariance  Cx  -*  oo  corresponding  to  no  a  priori  knowledge, 
we  obtain  the  simple  result  - 

P  =  Pmap  =  &2FTF 


(17) 


A3:  PROOF  THAT  GAVRIELIDES  FUNCTIONS  ARE 
THE  MINIMUM  NORM  VECTOR  POLYNOMIALS 

Consider  a  set  of  vector  functions  Fk(x)  obeying  - 


V  ■  Fk(x)  =  Zk(x) 

where  Zk{x )  are  the  wavefront  basis  functions  (Zernikes) 
and  which  satisfy  the  boundary  condition  - 


Fk(x)  -dl  =  0 


(41) 


(42) 


We  have  shown  in  section  3  of  this  report  that  the  modal  coefficients  of  the  wavefront 
expansion  can  then  be  estimated  by  direct  integration  as  - 


ak 


=  -  /  Fk{x) 
Jd 


V<f>(x)d2x 


(A3) 


Consider  adding  any  function  Ak{x)  to  Fk(x)  with  satisfying  the  same  boundary 

condition  A2  but  having  zero  divergence  (i.e.  V  •  Ak(x)  —  0).  It  is  then  clear  that  the 
functions  Fk(x)  are  not  unique  since  Fk(x)  +  Ak(x)  is  also  a  valid  solution. 

Gavrielides’  functions  which  we  will  denote  as  Gk(x)  are  chosen  to  be  equal  to  the 
gradient  of  an  associated  scalar  function  C- ( x_ ) •  Thus,  we  have  - 

Gk(x)  =  Vdk(x)  (M) 


Substitution  of  this  relationship  into  eqs  A1  and  A2.  yields  - 

V26(^)  =  Zk(x)  (A5) 

V(k(x)-dl  =  0  (-46) 

which  is  Poisson’s  equation  with  the  corresponding  boundary  condition.  Gavrielides  func¬ 
tions  are  obtained  through  solution  of  eq  A5  subject  to  the  boundaiy  condition  A6. 

Consider  a  set  of  vector  polynomials  obeying  eqs  A1  and  A2  given  by  Fk(x)  =  Gk(x)  + 
Ak(x)  where  the  Gk(x)  are  Gavrielides  functions  and  Ak\x)  are  an  arbitrary  set  of  functions 
obeying  eqs  A1  and  A2.  It  is  simple  to  show  that  the  noise  propagation  of  this  set  is 
governed  by  the  vector  norm  of  the  functions  - 

/  Fk(x)  ■  Fk(x)d2x  (-47) 

Jd 

The  noise  propagation  due  to  Gavrielides  functions  is  - 

Kg  =  I  Gk(x)  •  Gk(x)d2x 
Jd 


(-48) 


Using  eq.  A5,  standard  vector  identities  and  the  divergence  theorem,  this  may  be  expressed 


J  Zk(x)V2Zk(x)d2x  +  jZk(x)V£k{x)  -dl 

(A9) 

Substitution  of  eqs  A4  and  A5  into  A9  yield  - 

Ng  =  f  tk{x)Zk{x)d2x 

Jd 

(A10) 

The  noise  propagation  of  the  functions  Fk(x).  may  be  expressed  as 

- 

Nf  =  No  +  Na  +  2  /  Gk(x)  ■  Ak(x)d2x 

Jd 

(All) 

where  Na  is  given  by  - 

Na  =  /  At(x)  ■  Ak(x)d2x 

Jd 

(A12) 

The  integral  fD  Gk{x)  •  Ak{x)d2x  may  be  written  as  - 

Na  =  [  V^(x)  •  Ak(x)d2x 

Jd 

(A13) 

which  may  in  turn  be  expressed  as  - 

Na  =  J  V  •  Uk(x)Mx))d2x  -  £k(x)V  •  Ak(x)d2x 

(A14) 

The  second  integral  in  eq.  A14  is  identically  zero  since  the  functions  Ak(± 
0.  We  may  apply  the  divergence  theorem  to  the  first  term  m  eq.  A14  - 

■)  obey  V-Ajfe(x)  = 

J  V  •  (£k(x)Ak(x))d2x  =  j ik(x)Ak{x)  ■  dl 

(A15) 

but  since  the  Ak(x)  obey  the  boundary  condition  Ak(x)  ■  dl  -  0.  the  integral  fD  Gk(x)  ■ 

Ak(x)d2x  =  0  and  we  have 

Nf  —  Ng  +  N  a 

and  since  NA  is  by  definition  guaranteed  to  be  positive,  we  have 

(A16) 

Nf  Ng 

(A17) 

for  arbitrarv  Ak(x).  The  Gavriehdes  functions  are  thus  minimum  norm. 

Q.E.D. 

